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Abstract 

We present two new accurate and efficient method to compute the formal solution of the polarized 
radiative transfer equation. In this work, the source function and the absorption matrix are approxi- 
mated using quadratic and cubic Bezier spline interpolants. These schemes provide 2"^ and 3^^ order 
approximation respectively and don't suffer from erratic behavior of the polynomial approximation 
(overshooting). The accuracy and the convergence of the new method are studied along with other 
popular solutions of the radiative transfer equation, using stellar atmospheres with strong gradients 
in the line-of-sight velocity and in the magnetic-field vector. 

Subject headings: Line: profiles — Magnetic fields — Polarization — Radiative transfer — Stars: 
atmospheres 
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1. INTRODUCTION 

During the past decade, large scale numerical mod- 
els and inverse problem applications became powerful 
tools for diagnosing and understanding spectroscopic 
and spectropolarimetric observations. Massive applica- 
tions like 3D hydrodynamic (HD) simulations. Magnetic 
Doppler Imaging and data inversion of solar surface lay- 
ers stimulated interest in developing fast and robust for- 
mal solvers for the radiative transfer equation (RTE). 
These computationally demanding problems pose spe- 
cial requirements for the RTE solver: a large number 
of integrations needed to compute synthetic profiles per 
single iterations (time step or model adjustment). Fast 
convergence become particularly important as it allows 
achieving acceptable accuracy for the radiation field us- 
ing the same geometrical grid as used for hydrodynamics 
or inversion. 

The properties of RTE solver become even more critical 
for non-equilibrium modeling. The assumption of local 
thermodynamic equilibrium (LTE) makes level popula- 
tion densities defined by the local temperature of the 
model and decoupled from the radiation field. Only one 
integration of the RTE is needed to compute the emerg- 
ing intensity at each wavelength. However, in the non- 
LTE case (NLTE), evaluation of population densities is 
calculated using a self-consistent iterative method, nor- 
mally requiring two inte grations of the RTE per iteration 
(|01son fc Kunaszlll987[ ). Inaccuracies in the integration 
of the RTE on a coarser grid can lead to a slower con- 
vergence of the NLTE problem in the spectral synthesis 
and generally to a slower convergence of the inversion 
(in LTE and NLTE). The situation is worse when the 
magnetic field becomes important. 

The monochromatic radiative transfer equation for po- 
larized light can be expressed as: 



as 



(1) 



where / = (/, Q, /7, V)^ is the Stokes parameter vector, 
J= {ji^jQ^ju^jv)^ is the total emission vector, and K 
is the total absorption matrix. This matrix has seven in- 
dependent terms: 77/ is the total absorption of radiation; 



77Q, r]u and rjv describe the coupling of the intensity / 
with Q, U and V; and pg, pu and pv are cross-talk terms 
between Q, U and V due to mag neto-optical effects (see 
iLandi Degrinnocenti Landolfil [2QQ1 ) : 
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A number of advanced integration schemes, suitable 
for HD simulations and NLTE inversion s have bee n 
implemented for the unpolarized light (|Auerl l2QQ3f ). 
however, for polarized light it is commo n to use a 
short characteristics schemes assur ning linear (|Rees et all 
19891) or parabolic de pendence (|Trujillo Buenol I2QQ3I : 



Sampoorna et al]l2QQ8f ) of the source function with opti- 



c al path. 

iBellot Rubio et aH () 19981 ) proposed an efficient third- 
order method (LBR hereafter) that provides accurate 
results on coarse grids of depth-points for polarized 
light. This meth od has been used ex tensivel y to compute 
LTE inversions (jBellot Rubio et aI]l2 QQQ: S ocas-Navarro 
2011 ) and NLTE inversions (jde la Cruz Rodriguez et al 



20121 ). 

In this paper, we propose a new method to accurately 
integrate the polarized RTE, using Bezier interpolants. 
The paper is arranged as follows: an introduction to the 
DELO method and Bezier Splines are provided in §[2] and 
[3] respectively. For completeness, we introduce the solu- 
tion to the unpolarized RTE, using Bezier interpolants, 
in § m The quadratic and cubic Bezier solutions to the 
polarized RTE are presented, for the first time, in § [5l 
Our proposed solutions are tested against other meth- 
ods commonly used in radiative transfer studies in § [6l 
Finally, our conclusions are summarized in § [71 

2. THE DELO METHOD 

iRees et all ([19891 ) rewrite the polarized RTE using; the 
modified source vector {S = j/rji)^ dividing all terms in 
Equation ([T|) by the absorption coefficient 77/. This sim- 
ple change leaves the K matrix with all the diagonal 



elements normalized to one so the method got a name 
Diagonal Element Lambda Operator (DELO). 
Defining, 
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where 1 is the 4x4 identity matrix. The radiative trans- 
fer equation becomes: 



dl 



(4) 



with_the optical path defined as dr^ = rjids and 5^ = 
S-KL 

In our discrete grid of n depth points {k = 1, 2, 3, n), 
the solution to Equation (|4]) on the interval (r/e,r/e+i) is: 

r'T k+i 



I{Tk) = /(r/e+i)e" 



-5k 



<^-^^^5^{T)dT, (5) 



where 5k = drk. 

To integrate analytically Equation (j5j), iRees et aTl 
(|1989f ) assume that the equivalent source vector 5^ is 
linear with optical depth (conventionally referred to as 
DELO-linear). A quadratic dependence of 5^ with op- 
tical depth would ideally allow to improve the conver- 
gence, however, it becomes unstable in the presence of 
non- linear gradients (Murphy 1990). 

3. THE BEZIER INTERPOLANTS 

lAuerl ()2QQ3f ) provides a thorough summary of advanced 
interpolants that could be used to integrate the (implic- 
itly assumed) unpolarized radiative transfer equation. 

3.1. Bezier quadratic interpolant 
Defining normalized abscissa units in the interval 

(Xfe,X/e+i), 

(6) 



the quadratic Bezier interpolant can be expressed as: 

f{x) = (1 - ufvk + Vk+iu^ + 2u{l -u)-C, (7) 

where yk and yk-\-i represent the node values of the func- 
tion that is being interpolated and C is the Bezier control 
point. The latter can be expressed using the derivative 
at Xk or x/c+i. 



C^{xk) = yk + Y^'(xfe), 
C^{xk+i)=yk+i - ^y\xk+i), 



(8) 
(9) 



with hk = Xk-\-i — If both and can be com- 
puted, it is desirable to take the mean: 



C= (CVC^)/2, 

resulting in a more "symmetric" spline. 
Defining, 
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Fig. 1. — Comparison of 3 interpolation schemes. The black 
circles indicate the know values of the function that is interpolated 
using quadratic Bezier spline, cubic Bezier spline and piece-wise 
parabolic fit. 



accurate numerical derivatives are calculated at the 
node points {xi^2,...,n) with a scheme that sup- 
presses overshooting beyond the node function values 
(Frit sch fc Butland.1984) . The Bezier derivative is given 
by: 

y'i-^) = . . . — , (14) 



: • d, 



+ 1/2 



+ (1 - a) 'di_i/2' 



if dk-i/2 ' dk-\-i/2 > and it is set to otherwise. 

3.2. Bezier cubic interpolant 

Alternatively, the cubic Bezier interpolant can be ex- 
pressed as: 

f{x) =(1 - ufyk + yk+iu^^ 

3u{l - ■ E ^ 3u^{l - u) ■ F, 



(15) 



where E and F are control points, similar to those de- 
fined in Equations (|8]) and (|9]), but each of them is placed 
in a different location: 



E{xk) = yk + Y^'(^fe), 

F(Xfe+l)=7//e+l - ■y?/'(x/e+l). 



(16) 
(17) 



Fig. [T] illustrates the behavior of three interpolation 
schemes, when the function has fast changing gradient. 
Here, we compare a parabolic polynomial, the quadratic 
Bezier splines (Bezier 2) and the cubic Bezier splines 
(Bezier 3) described in this section. 

The behavior of the parabolic fit is erratic in the vicin- 
ity of intervals with large change of gradient, showing 
dramatic overshooting peaks. The Bezier splines consid- 
ered here stays between t he pai r of bo undary data values 
in each interval. If fact, lAuerl (|2QQ3f ) already suggested 
that the Bezier splines (and the Hermitian interpolant) 
are particularly suitable for solving the radiative trans- 
fer equation in form ^ as they do not produce spurious 
maxima and minima. 

From x = — ltox = l the curve increases smoothly, 
and the three interpolation schemes considered here pro- 
duce very similar values. 
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4. BEZIER INTEGRATION OF THE SCALAR RTE 

In this section, we derive the formal solution_of the 
radiative transfer equation for unpolarized hght (K = 0). 
In this case Equation ^ becomes a scalar equation: 



dJ 

dr 



(18) 



where / is the intensity and S is the unpolarized source 
function. Equation (j5|) is easily integrated by approxi- 
mating S with any of the Bezier interpolants described 
in Sec. [H 

4.1. Quadratic Bezier integration 

If the quadratic Bezier interpolant is used to approx- 
imate the source function, the solution to Equation (|5]) 
is 



/(r/e) = /(r/e+i) 



+ akSk + PkSk+i + IkCk, (19) 



where Pk and are coefficients from the integral 
that only depend on Sk = r/c+i — r^: 
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To calculate 6^, the opacity can be integrated ana- 
lytically by using Bezier splines to approximate r]i as 
a function of depth. Further details can be found in Ap- 
pendix [Al 

Note, that for small Sk it is wise to use Taylor expan- 
sion for the exponents in the right-hand-side to avoid di- 
vision of vanishingly small quantities (see Appendix [B]) . 

Equation (p!9|) has the form Ik = A - Ik-\-i -\- B. In the 
case of a stellar atmosphere we can sequentially evaluate 
outgoing intensity at any depth point starting from the 
boundary condition In = Sn at the bottom. The control 
point Ck is computed using Equation (p!Q|) for all inner 
points of the integration domain. For the top point only 
one approximation for the control point (C^) is available. 

4.2. Cubic Bezier integration 

If the cubic Bezier interpolant is used to approximate 
the source function, the solution to Equation ([5]) is 

-Sk 

, ^ . „ (20) 
7kEk 



I{Tk)=I{Tk+i)-e-'"'+akSk 



o^k, $k, 7fe, Cfc are coefficients from the integral that 
only depend on 4: 
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The same recommendations apply here for small Sk- 

5. BEZIER INTEGRATION OF THE VECTOR RTE 

5.1. Quadratic Bezier integration 

For the polarized case (K 7^ 0), the solution to Equa- 
tion (j5j) can formally be presented in the form similar to 
Equation ([T9)) : 



/(rfc) = /(r/e+i) • e ' 

^Pk5^{Tk+l) 



+ c^fe^(r/c) 
• IkCk. 



(21) 



In Sec. m Equation ([H]) we show how to compute the 
derivatives of S needed for evaluating the control point 
C. In the polarized case, however, derivatives of J, S 
and K must be computed. The presence of Ik in the de- 
nominators of both and (through Equation (fT4|) ) 
introduces non-linearity that kills a simple recurrence re- 
lation available in the scalar case: Ik = A ■ Ik+i + B. 

An elegant solution to this problem is actually provided 
by the radiative transfer equation: 



(22) 

— — *^/c+l ^ (23) 

Using these expressions, the control points can be re- 
written as: 



■-Ik — Sk ^ Kfe/fe, 

'- Ik-\-l — Sk-\-l +K/e+i//c+i. 
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and: 
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(25) 



Note that Cj. and are vectors, whereas is a matrix. 

The derivatives of the modified opacity matrix and 
the source vector must be computed according to the 
recipe for monotonicity given by Equations (pT]) - (fT4|) to 
prevent the overshooting. 

This convenient splitting of permits re-writing 
Equation ([2T]) for unknown Stokes vector in point /c as a 
system of four linear equations {A-Ik = B): 



[t^ak'Kk 



Ik. 



with, 



^k 
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akSk + 



(26) 



(27) 



Note, that the matrix and the right-hand-side in Equa- 
tion [26l include only known quantities of absorption ma- 
trix, source and Stokes vectors. Note also, that solving 
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the system of linear equations (|26]) directly instead of 
trying to invert A is both faster and more robust. 

5.2. Cubic Bezier integration 

Similarly to the quadratic Bezier integration, the solu- 
tion to Equation (j5j) is formally similar to Equation (|2Q|) : 

+ /3fc5^(rfc+i) + %Ek + Ckh- (28) 

Given the formal similarity of this equation with the 
quadratic case, where also two control points are used, 
the algebra needed to write the problem as a linear sys- 
tem of equations (Equation (|26|) ) is almost identical. 



-K/e ) Sk 



(KfeKfc + K'fe + Kfe) + 



_ (29) 



and: 



1 T^K/e+i I Sk-\-l 



4, 
3 



/c+l" 



+ 



— (Kfc+iKfc+i +K5j_,_i +Kfc+i) 



(30) 



n dSHl) can 
; I{rk) on 

(l + akKk - %ek) Ik = 



Equation (|28|) can be re-arranged to group all the terms 
containing I{rk) on the left-hand side of the equation. 



-5k 
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(31) 



Again we advise to solve the linear system of equations 
in Equation (|3T]) . instead of inverting the A matrix. 

6. NUMERICAL CALCULATIONS 

We use a modified version of the radiative transfer 
code Nicole (Socas-Navarro et al. in prep), to test 
the numerical accuracy of the DELO-Bezier methods 
for polarized light. A snapshot from a 3D MHD nu- 
merical simulation is used to carry-out the calculation 
of synthetic full- stokes profiles. It covers a physical 
range of 16.6 x 8.3 x 15.5 Mm, extending from the up- 
per convection zone to the lower corona (from 1.5 Mm 
below to 14 Mm above average optical depth unity at 
5000 A). The simulation has an average magnetic field 
strength of 150 G in the phot osphere. This p articu - 
lar snapshot has been used by iLeenaarts et aTl ()2009l ): 
Ide la Cruz Rodriguez et al ] (HoH) so further details on 
simulations can be found there. 

6.1. The NLTE problem 

As described in Ide la Cruz Rodriguez et all (120121 ). the 
Nicole code solves the NLTE problem for unpolarized 



light (| Socas-Navarro fc Trujillo Buenol 119971 ) and, once 
the atom population densities are known, a polarized 
formal solution is calculated (the polarization-free ap- 
proxi mation, see iTrujillo Bueno fc Landi Degrinnocentil 
Il996f ). The Ca II atom model used in this work consists 
of five bound levels plus a continuum. 

Additionally, the velocity- free approximation, pre- 
vious ly used to comput e the data inversions in 
Socas-Navarro et al.l ()2000f ): Ide la Cruz Rodriguez et aTl 
(j2012[ ). is not used in this study. Therefore the popula- 
tion densities are fully consistent with the strong velocity 
gradients present in our models. 

For consistency, we have implemented an accelerated 
local lambda operator based in the unpolarized solution 
described in Sec. IH 

The process of constructing the approximate lambda 
operator at any given point /c, can be idealized by set- 
ting the source function to unity in that point, and 
zero otherwise. The operator is constructed using the 
SI that remain after this opera- 




According to Equation (p!Q|) . the average control point 
C has an explicit dependency on C^{Sk) and Cl{Sk-\-i)- 
Using the average of the two Cs improves stability and 
accuracy of the Bezier interpolant, although the lambda 
operator becomes intrinsically non local. 

A simple recipe for reducing the number of itera- 
tions for the NLTE problem is to make the approximate 
lambda operator strictly local by setting C = C^. Rare 
cases of minor overshooting are suppressed by forcing 
the control point mmjSk, Sk-\--[) < Ck < mdix(Sk,S k-\-i)- 
Similar strategies are found in iHavek et al.l (|2010f ) and 
iHolzreuter fc Solankil l\2012\ ). 

We used the approximate lambda operator for solving 
the NLTE problem in the form: 



(32) 



where r^^^ is the r atio between the lin e opacity and the 
total opacity (see iRvbicki fc Hummerl IT99 1. for further 
details). 

Fig. [2] illustrates the full-Stokes images at the surface 
of a 8.3 X 8.3 Mm patch computed at —10 mA from the 
core of the Ca II A8542 line. Columns of panels (from left 
to right) show the monochromatic images obtained using 
DELO-parabolic, DELO-linear, quadratic DELO-Bezier 
and cubic DELO-Bezier, respectively. We used identical 
NLTE level populations to compute the polarized for- 
mal solution with each algorithm. Thus, the differences 
between panels in different columns reflect the numeri- 
cal properties (convergence and stability) of the methods 
compared. 

The DELO-parabolic method seems to work in most of 
the pixels, but it fails notoriously to produce an accurate 
solution at certain areas which have been indicated on 
the panels using red markers. The artifacts here can be 
quite large, and some pixels show even negative values 
in Stokes /. 

The maximum difference in Stokes / between DELO- 
linear and DELO-Beziers at this wavelength is less than 
1% of the continuum intensity. This is an expected result, 
given that the vertical grid spacing of 3D snapshots is 
so fine that even a linear scheme produces an accurate 
solution. 
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Fig. 2. — Monochromatic full-Stokes images computed at the core of the Ca II A8542 line. From left to right, the images are computed 
using the DELO-parabolic, DELO-linear, DELO-Bezier and cubic DELO-Bezier, respectively. Stokes /, Q, U and V are represented from 
top to bottom respectively. Regions with artifacts are marked in the panels computed using the DELO-par method using small red squares. 



This example only shows the stability of the higher- 
order Bezier methods, but it hardly shows any advantage 
over a linear scheme. 

In Sec. 16.21 we test the converges and stability of these 
methods using ID models where strong gradients in the 
magnetic-field and line-of-sight (l.o.s) velocity are intro- 
duced. 

6.2. Numerical accuracy 

To assess the accuracy of the formal solvers, we have 
created four atmospheric models using a fine grid of 
depth-points with 198 points-per-decade, equidistantly 



placed from logrsoo = —6.9 to logrsoo = 2. 

The four atmospheric models share the same temper- 
ature, electron pr essure, and micro-tur bulence from the 
VALC-III model (jVernazza et al.l[T98TI ). 

We have created a complicated ad-hoc magnetic field 
vector and l.o.s velocity component, which are illustrated 
in Fig. [3l along with the VALC-III temperature. To cre- 
ate a demanding topology, the horizontal component of 
the magnetic-field follows a spiral centered on the line- 
of-sight. 

Considering the quantities described above, our set of 
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Fig. 3. — Physical quantities from our prescribed ID stellar 
model. The top panel shows the temperature stratification as a 
function of optical depth. The middle panel shows an ad-hoc pro- 
file used to define vi_o.s our calculations. The bottom panel il- 
lustrates the three components of the magnetic field. To illustrate 
the effects of having a poor sampling in the vertical dimension, we 
have marked in the first two panels, the equivalent model assuming 
1 point-per-decade (blue circles), 3 points-per-decade (red crosses) 
and 198 points-per-decade (black dots). Note that the z-axis is 
parallel to the line-of-sight. 

four models is summarized as follows: 

1. The first case is a model with constant l.o.s veloc- 
ity {vi,o.s =0 km s~^) and constant magnetic-field 
B = (600, -700, 800) G for ah depth-points. 

2. In the second case, the constant line-of-sight veloc- 
ity is replaced with our ad-hoc profile. 

3. The third case has constant velocity but the com- 
plicated magnetic-field vector described in Fig. [3l 

4. Complex velocity and magnetic-field profiles are 
used simultaneously in the fourth case. 

For each of these four cases, we have computed the 
absorption matrix and the source function in the orig- 
inal grid of depth-points for the Ca II A8542 and the 
Fe I A6302 lines. These two lines have been extensively 
used as atmospheric diagnostics in solar applications. 
The former is sensitive to a vast range of height includ- 
ing photospheric and chromospheric response, although 



it has an relatively low Lande factor (^eff = l-l)- The lat- 
ter is formed in the solar photosphere, but it has a much 
higher sensitivity to the magnetic-field (^eff = 2.5). 

Our test consists of computing the formal solution to 
the polarized RTE using a pre-computed absorption ma- 
trix and source function. To study the converge prop- 
erties of the integrators we simply drop intermediate 
depth-points from the original grid. Thus we separate 
the level population calculations from the RTE solution. 
We also noticed that our level population calculations for 
Ca II start deviating noticeably from the correct solution 
when the number of depth-points is below 10 points-per- 
decade. 

As a reference here we use the emerging Stokes vector 
computed on the original grid of 198 points-per-decade, 
using an adaptive fourth-o rder Runge-Kutta solver (see 
iLandi Degrinnocentil 119761 ). which ensures that the pre- 
cision of 10~^^ as achieved at every step. 

Fig. m shows a set of Stokes /, Q, /7 and V pro- 
files computed for the complex of our ad-hoc models 
(model 4) using a 2 points-per-decade grid. The reference 
grey profile has been computed using the Runge-Kutta 
method and the original dense grid of depth-points. The 
main discrepancies occur at the chromospheric NLTE 
core of the Ca II A8542 line. At those wavelengths, 
the monochromatic depth-scale presents large irregular 
jumps in optical-depth. 

For example, in Stokes / the LBR and DELO-linear 
solutions present an excessively strong line core, whereas 
the DELO-parabolic method produces a line core in emis- 
sion. The DELO-Bezier profiles tightly trace the ref- 
erence profile at all wavelengths, and only a noticeable 
small deviation is present at the very inner core of the 
line. 

The results of our calculations are summarized in Fig. [5] 
and Fig. El for the Ca II A8542 line and the Fe I A6302 
line respectively. The line profiles are compared with the 
reference case, and the largest discrepancy over the entire 
profile (in absolute value) is plotted for each integration 
method as a function of the number of depth-points per 
decade. 

This error measurement is sensitive to outliers at a 
single wavelength. Therefore, the curves in Fig. I5] and 
Fig. ini may not be monotonous functions of the number 
of grid points. Still, we prefer this metric as it gives 
a robust estimate of accuracy, convergence speed and 
stability for our solvers. 

Ideally, we would expect that any discrepancies shown 
by all the methods considered here, would disappear 
when the density of depth-points is large enough (all 
methods should converge to the same accurate solution). 
The errors are expected to be large when the density 
of depth-points is low, especially for the highest order 
methods, but these should achieve a more accurate solu- 
tion than the lower order methods when the density of 
depth-points is relatively high. 

In Fig. 13 our results show that at the absolute mini- 
mum of 1 points-per-decade, it is hard to beat the DELO- 
linear solution, which keep the errors under control de- 
spite the strong gradients and poorly-sampled stellar at- 
mosphere. However, the DELO-Bezier solutions provide 
almost as accurate results except in Stokes Q where the 
error is half an order of magnitude higher. 
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Fig. 4. — Emerging intensity vector at the Ca II A8542 line. The 
profiles are computed using our fourth model which contains gra- 
dients in the magnetic-field and line-of-sight-velocity. The model 
here has a vertical resolution of of 2 points-per-decade. The ref- 
erence profile is computed using a Runge-Kutta solver using the 
original grid of 198 points-per-decade. For readability, only the 
cubic DELO-Bezier solution is plotted. 

The situation changes drasticahy, as a few more depth- 
points are included in the calculations. Between five 
and thirty points-per-decade, the DELO-Bezier solutions 
clearly out-perform all the other solutions for the four 
scenarios that we have prepared. Most radiative transfer 
computations are carried out within this range of grid 
densities. This also matches hydrodynamic grid density 
typically used for environments where radiation carries 
important fraction of energy. 

The DELO-parabolic needs a large number of depth- 



points to match the performance of all the other meth- 
ods, as was expected. 

The calculations for the Fe I A6302 line, show a 
smoother behavior than in the previous case and all cases 
seem to reach a stable level of accuracy at approximately 
30 points-per-decade. The inclusion of more points above 
this value, seems to increase the accuracy marginally and 
very inefficiently. 

The DELO-Bezier demonstrated fastest convergence 
reaching below 10 % for the least dense grids in Stokes /. 
This conclusion is true for all our models. The accuracy 
of all methods is comparable as they all reach the same 
level of precision above 30 points-per-decade. 

For both spectral lines, the performance of the LBR 
method is close to both the DELO-linear and DELO- 
parabolic methods, and it only seems to be slightly less 
accurate when strong gradients in the line-of-sight veloc- 
ity are present. 

The accuracy achieved by the DELO-Bezier methods 
is almost identical. The quadratic solution marginally 
outperforms the cubic one when the number of depth- 
points per decade is 1 or 2, whereas the cubic version is 
slightly better at the Ca II A8542 line, within the range 
2 < Ndepth < 10 points-per-decade. 

Computationally, the LBR and DELO-Bezier solutions 
are similarly demanding given that accurate derivatives 
of the source function and of the absorption matrix are 
required. In comparison, the DELO-linear and DELO- 
parabolic solutions are almost two times faster, because 
no derivatives are required. 

Note also that, in the first and last intervals of the grid, 
only one control-point can be computed. Therefore, es- 
pecial care must be taken when the cubic method is used. 
The simplest solution is to use the quadratic method in 
those two intervals, given that only one control point is 
strictly necessary. 

7. DISCUSSION AND CONCLUSIONS 

In this work we present a new scheme to integrate the 
polarized radiative transfer equation using Bezier splines. 
The solvers are a general ization of the advanced strategy 
proposed by lAuerl (|2003l ) for the unpolarized case. Our 
second-order and third-order integration methods pre- 
serve the stability of the DELO-linear solution with the 
benefits of faster convergence and higher order accuracy. 

The advantages of the new integration scheme are il- 
lustrated by comparison with other methods commonly 
used to solve the polarized RTE. We show that the new 
DELO-Bezier schemes outperform all other algorithms 
when the density of depth-points is low and the struc- 
ture of the medium includes steep gradients. 

The new methods will be beneficial for various applica- 
tion from reconstructing atmospheric structure (s) using 
observational data (data inversion in solar physics. Mag- 
netic Doppler Imaging of stars) to radiative magneto- 
hydrodynamics. E.g., in case of solar data inversion the 
new solvers will allow using less denser depth grid im- 
proving the stability of inversion, which is particularly 
important for strong lines, like the Ca II H,K and the 
IR triplet. Given the vast formation range of these lines 
(larger than 1000 km), it is hard to find an optimal grid 
of depth points at all frequencies in the line. 

The Institute of Theoretical Astrophysics (Univer- 
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Fig. 5. — Maximum error for the Ca II A8542 line as a function of the number of points-per-decade (at all wavelengths). From top to 
bottom, the panels show the errors for Stokes /, Q, U and V, respectively. From left to right, each column corresponds to each of the four 
models used in our calculations. Each colored line corresponds to one of the formal solutions used in our calculations: LBR (solid purple), 
DELO-linear (dotted blue), DELO-parabolic (dashed orange), quadratic DELO-Bezier (solid green) and cubic DELO-Bezier (dashed black). 
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Fig. 6. — Maximum error for the Fe I A6302 line as a function of the number of points-per-decade (at all wavelengths). From top to 
bottom, the panels show the errors for Stokes /, Q, U and V, respectively. From left to right, each column corresponds to each of the four 
models used in our calculations. Each colored line corresponds to one of the formal solutions used in our calculations: LBR (solid purple), 
DELO-linear (dotted blue), DELO-parabolic (dashed orange), quadratic DELO-Bezier (solid green) and cubic DELO-Bezier (dashed black). 
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APPENDIX 
A. BEZIER INTEGRAL OF THE OPACITY 

In many radiative transfer applications, a conversion of the depth-scale is desirable (i.e, from z to logrsoo). To 
carry out this conversion, the opacity = rj) is integrated along the depth-scale. When the grid of depth-points 
is sufficiently dense, a trapezoidal integration would suffice. This wide-spread method becomes inaccurate for coarse 
depth-scale grids. 

We advice to approximate the opacity with a quadratic Bezier spline (Eq.[7j), and integrate analytically this function. 
Once again, the integral can be expressed as a set of coefficients that multiply the values of the opacity that define the 
integration interval: 

drk{iy) = r/e+i - =(x/e+i - Xk) / [vki^^) * (1 - + Vk+i{^) * + 2u{l - u)C] du = 

•^0 (Al) 

B. TAYLOR EXPANSION OF THE BEZIER INTEGRAL COEFFICIENTS 
In §|4]we advise to use a Taylor expansion of the exponential term in Equation (j5j), when 8^ is small: 

e ^^^^1-4^/+^ 

The resulting integration coefficients for the quadratic Bezier interpolant are: 
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The integration coefficients for the cubic Bezier interpolant are: 
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